% Forecast Results

% Descriptive Statists
clear all;
small = 1.0e-10;

% Model Name
mod_name = '_baseline';
% mod_name = '_fgamma05';
% mod_name = '_fgamma15';
% mod_name = '_q_8';
% mod_name = '_q_20';

% -- File Directories   
figdir = 'fig/';
outdir = 'out/';
matdir = 'mat/';
mcmcdir = 'mcmc/';
figdir_fcst = 'fig_fcst/';

fstr = [matdir 'country_code']; load(fstr);
fstr = [matdir 'country']; load(fstr);
n_c = size(country_code,1);

fstr = [matdir 'yfit_mat'];load(fstr);  %Trends .. for plotting

% load Paths
fstr = [mcmcdir 'path_F_draws' mod_name];load(fstr);
fstr = [mcmcdir 'paths_U_draws' mod_name];load(fstr);
fstr = [matdir 'yp']; load(fstr);
fstr = [matdir 'pop']; load(fstr);

% Rename;
u_draws = paths_U_draws;
f_draws = path_F_draws;
clear paths_U_draws;
clear path_F_draws;

T = 118;  % Number of insample obs

% Compute Growth Rate Paths
U_gr_paths_50 = squeeze(u_draws(T+50,:,:)-u_draws(T,:,:))/50;
U_gr_paths_100 = squeeze(u_draws(T+100,:,:)-u_draws(T,:,:))/100;
f_gr_paths_50 = (f_draws(T+50,:)-f_draws(T,:))/50;
f_gr_paths_100 = (f_draws(T+100,:)-f_draws(T,:))/100;
Y_gr_paths_50 = U_gr_paths_50 + repmat(f_gr_paths_50,n_c,1);
Y_gr_paths_100 = U_gr_paths_100 + repmat(f_gr_paths_100,n_c,1);

% Construct Growth Rate Prediction intervals for each country
% Save Results 
outfile_name = [outdir 'Table_3' mod_name '.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,'Forecast Growth Rates (PAAR) \n');
fprintf(fileID,',,50Years,,,,,,,100Years \n');
pct = [0.05 1/6 0.50 5/6 0.95]';
str1 = ['Mean,.05,.16,.50,.84,.95'];
fprintf(fileID,[',,' str1 ',,' str1 '\n']);
fprintf(fileID,'F,,');
tmp = pctile_mean(f_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(f_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f','\n');

% Use population weights for growth rates
fprintf(fileID,'All,,');
pop_g = pop(end,:)';     % All countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50;
g_gr_paths_100 = wght'*Y_gr_paths_100;
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

% Load Lists
load([matdir 'oecd_list']);
i_oecd = zeros(n_c,1);
nn = size(oecd_list,1);
for j = 1:nn;
    ss = char(oecd_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_oecd(jj) = 1;
    end;
end;

ii = i_oecd==1;
fprintf(fileID,'\n OECD,,');
pop_g = pop(end,ii==1)';     % OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_oecd==0;
fprintf(fileID,'\n Non-OECD,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

% Load WEO lists of countries
fstr1 = ['List_of_countries_weo.xlsx'];
weo_lists_ii = readmatrix(fstr1,'Range','C3:L115');
weo_country_name = readmatrix(fstr1,'Range','B3:B115','OutputType','char');
weo_country_code = readmatrix(fstr1,'Range','A3:A115','OutputType','char');

list_cc_adv_econ = weo_country_code(weo_lists_ii(:,1) == 1);
list_cc_ea = weo_country_code(weo_lists_ii(:,2) == 1);
list_cc_g7 = weo_country_code(weo_lists_ii(:,3) == 1);
list_cc_em_dev = weo_country_code(weo_lists_ii(:,4) == 1);
list_cc_asean5 = weo_country_code(weo_lists_ii(:,5) == 1);
list_cc_em_dev_asia = weo_country_code(weo_lists_ii(:,6) == 1);
list_cc_em_dev_europe = weo_country_code(weo_lists_ii(:,7) == 1);
list_cc_latin = weo_country_code(weo_lists_ii(:,8) == 1);
list_cc_me = weo_country_code(weo_lists_ii(:,9) == 1);
list_cc_ss_africa = weo_country_code(weo_lists_ii(:,10) == 1);

g_list = list_cc_adv_econ;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_adv_econ = i_g;

g_list = list_cc_ea;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_ea = i_g;

g_list = list_cc_g7;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_g7 = i_g;

g_list = list_cc_em_dev;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_em_dev = i_g;

g_list = list_cc_asean5;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_asean5 = i_g;

g_list = list_cc_em_dev_asia;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_em_dev_asia = i_g;

g_list = list_cc_em_dev_europe;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_em_dev_europe = i_g;


g_list = list_cc_latin;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_latin = i_g;


g_list = list_cc_me;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_me = i_g;

g_list = list_cc_ss_africa;
i_g = zeros(n_c,1);
nn = size(g_list,1);
for j = 1:nn;
    ss = char(g_list(j));
    jj = colnumber(ss,country_code);
    if jj > 0;
      i_g(jj) = 1;
    end;
end;
i_ss_africa = i_g;

ii = i_adv_econ==1;
sum(ii)
fprintf(fileID,'\n Advanced Economies,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_ea==1;
sum(ii)
fprintf(fileID,'\n Euro Area,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_g7==1;
sum(ii)
fprintf(fileID,'\n G7,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_em_dev==1;
sum(ii)
fprintf(fileID,'\n Emerging and Developing Economies,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_em_dev_asia==1;
sum(ii)
fprintf(fileID,'\n Emerging and Developing Asia),,')
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_asean5==1;
sum(ii)
fprintf(fileID,'\n ASEAN-5,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');


ii = i_em_dev_europe==1;
sum(ii)
fprintf(fileID,'\n Emerging and Developing Europe,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_latin==1;
sum(ii)
fprintf(fileID,'\n Latin America and the Caribbean,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
 
ii = i_me==1;
sum(ii)
fprintf(fileID,'\n Middle East and Central Asia ,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');

ii = i_ss_africa==1;
sum(ii)
fprintf(fileID,'\n Sub-Saharan Africa,,');
pop_g = pop(end,ii==1)';     % Non-OECD countries
wght = pop_g/sum(pop_g);
g_gr_paths_50 = wght'*Y_gr_paths_50(ii==1,:);
g_gr_paths_100 = wght'*Y_gr_paths_100(ii==1,:);
tmp = pctile_mean(g_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');
tmp = pctile_mean(g_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%5.2f',',,');


% Compute and Save Country-specific forecasts
% Construct Growth Rate Prediction intervals for each country
% Save Results 
outfile_name = [outdir 'Table_10_detail_byCountry' mod_name '.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,'Forecast Growth Rates (PAAR) \n');
fprintf(fileID,',,50Years,,,,,,,100Years \n');
pct = [0.05 1/6 0.50 5/6 0.95]';
str1 = ['Mean,.05,.16,.50,.84,.95'];
fprintf(fileID,[',,2017Value,,' str1 ',,' str1 '\n']);
fprintf(fileID,'F,,');
tmp = mean(f_draws(T,:)');
fprintf(fileID,'%7.4f ,,',tmp);
tmp = pctile_mean(f_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%7.4f',',,');
tmp = pctile_mean(f_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%7.4f','\n');

fprintf(fileID,'\n\n');
for ic = 1:n_c;
    str = char(country_code(ic));
    fprintf(fileID,[str ',,']);
    tmp1 = squeeze(u_draws(T,ic,:));
    tmp = mean(f_draws(T,:)') + mean(tmp1');
    fprintf(fileID,'%5.2f ,,',tmp);
    tmp = pctile_mean(Y_gr_paths_50(ic,:)',pct)';
    prtmat_comma(100*tmp,fileID,'%5.2f',',,');
    tmp = pctile_mean(Y_gr_paths_100(ic,:)',pct)';
    prtmat_comma(100*tmp,fileID,'%5.2f',',,');
    
    % Construct Median Log-level after 100 years
    ytmp = squeeze(u_draws(end,ic,:))+f_draws(end,:)';
    med_100 = median(ytmp);
    fprintf(fileID,'%5.2f \n',med_100);
end;
 


% Compute and Save Country-specific forecasts
% Construct Growth Rate Prediction intervals for each country
% Save Results 
outfile_name = [outdir 'Table_3_detail_byCountry' mod_name '.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,'Forecast Growth Rates (PAAR) \n');
fprintf(fileID,',,50Years,,,,,,,100Years \n');
pct = [0.05 1/6 0.50 5/6 0.95]';
str1 = ['Mean,.05,.16,.50,.84,.95'];
fprintf(fileID,[',,2017Value,,' str1 ',,' str1 '\n']);
fprintf(fileID,'F,,');
tmp = mean(f_draws(T,:)');
fprintf(fileID,'%7.4f ,,',tmp);
tmp = pctile_mean(f_gr_paths_50',pct)';
prtmat_comma(100*tmp,fileID,'%7.4f',',,');
tmp = pctile_mean(f_gr_paths_100',pct)';
prtmat_comma(100*tmp,fileID,'%7.4f','\n');

fprintf(fileID,'\n\n');
for ic = 1:n_c;
    str = char(country_code(ic));
    fprintf(fileID,[str ',,']);
    tmp1 = squeeze(u_draws(T,ic,:));
    tmp = mean(f_draws(T,:)') + mean(tmp1');
    fprintf(fileID,'%5.2f ,,',tmp);
    tmp = pctile_mean(Y_gr_paths_50(ic,:)',pct)';
    prtmat_comma(100*tmp,fileID,'%5.2f',',,');
    tmp = pctile_mean(Y_gr_paths_100(ic,:)',pct)';
    prtmat_comma(100*tmp,fileID,'%5.2f',',,');
    
    % Construct Median Log-level after 100 years
    ytmp = squeeze(u_draws(end,ic,:))+f_draws(end,:)';
    med_100 = median(ytmp);
    fprintf(fileID,'%5.2f \n',med_100);
end;


% Save Some Results -- used below
post_m_perc_save_multivariate = NaN(6,2,n_c);
for ic = 1:n_c;
    post_m_perc_save_multivariate(:,1,ic) = pctile_mean(Y_gr_paths_50(ic,:)',pct);
    post_m_perc_save_multivariate(:,2,ic) = pctile_mean(Y_gr_paths_100(ic,:)',pct);
end;
post_f_perc_save = NaN(6,2);
post_f_perc_save(:,1) = pctile_mean(f_gr_paths_50(1,:)',pct);
post_f_perc_save(:,2) = pctile_mean(f_gr_paths_100(1,:)',pct);


% Figure 5 -- Multivariate Prediction intervals by country
M50 = squeeze(post_m_perc_save_multivariate(2:6,1,:));
M100 = squeeze(post_m_perc_save_multivariate(2:6,2,:));

% Sort by average YP over last 10 years
yp_avg = mean(yp(end-4:end,:))';
[~,II] = sort(yp_avg);
country_code_srt = country_code(II);
M50 = M50(:,II)';
M100 = M100(:,II)';

% Error bands for plotting
trend = (1:1:n_c)';
trendp = trend+0.1;


% List country names for future use;
outfile_name = [outdir 'plot_univariat_multivariate_fcst_intervals_country_names' mod_name '.csv'];
fileID = fopen(outfile_name,'w');
for i = 1:n_c;
    fprintf(fileID,'%3i ,',i);
    str = char(country_code_srt(i));
    fprintf(fileID,[str '\n']);
end;

figure_5;


% Figure 6_a
% Read in Univariate results are compute predictive percentiles
U50_pct = NaN(n_c,5);
U100_pct = NaN(n_c,5);
% Read in UMs values
umdir = ['/Users/mwatson/Dropbox/Shared_Folders/SCC/output_matlab/ndata/'];
fstr = [umdir 'single_paths_baseline']; load(fstr);

% Compute univariate forecast intervals and save them
for ic = 1:n_c;
    pc = squeeze(paths(:,ic,:));
    g50 = 100*(pc(T+50,:)-pc(T,:))/50;
    tmp = pctile(g50',pct)';
    U50_pct(ic,:) = tmp;
    g100 = 100*(pc(T+100,:)-pc(T,:))/100;
    tmp = pctile(g100',pct)';
    U100_pct(ic,:) = tmp;
end;

% Save Univariate percentiles to a file
outfile_name = [outdir 'Univariate_Forecast_Intervals.csv'];
fileID = fopen(outfile_name,'w');
fprintf(fileID,'Country 50 the 100 year percentiles \n');
for ic = 1:n_c;
    str = char(country_code(ic));
    fprintf(fileID,[str ',,']);
    prtmat_comma(U50_pct(ic,:),fileID,'%5.3f',',,');
    prtmat_comma(U100_pct(ic,:),fileID,'%5.3f',',,');
    
    % Construct Median Log-level after 100 years
    ytmp = squeeze(paths(end,ic,:));
    med_100 = median(ytmp);
    fprintf(fileID,'%5.2f \n',med_100);
    
end;

M50 = squeeze(post_m_perc_save_multivariate(2:6,1,:));
M100 = squeeze(post_m_perc_save_multivariate(2:6,2,:));

% Sort by average YP over last 10 years
yp_avg = mean(yp(end-4:end,:))';
[~,II] = sort(yp_avg);
country_code_srt = country_code(II);
M50 = M50(:,II)';
M100 = M100(:,II)';
U50 = U50_pct(II,:)/100;
U100 = U100_pct(II,:)/100;

% Try error bands for plotting
trend = (1:1:n_c)';
trendp = trend+0.25;

figure_6_a;

%   Figure 4 .. individual countries
% Plot Individual figures for each country (save for appendix)
ffcst_pct = NaN(T+100,5);
for t = 1:T+100;
    ffcst_pct(t,:) = pctile(f_draws(t,:)',pct');
end;
yfcst_m_pct = NaN(T+100,n_c,5);
yfcst_u_pct = NaN(T+100,n_c,5);
for ic = 1:n_c;
    U_path = squeeze(u_draws(:,ic,:));
    Y_path = U_path + f_draws;
    Y_u_path = squeeze(paths(:,ic,:));
    for t = 1:T+100;
      yfcst_m_pct(t,ic,:) = pctile(Y_path(t,:)',pct');
      yfcst_u_pct(t,ic,:) = pctile(Y_u_path(t,:)',pct');
    end;
end;

calvec_p = (1900:2117)';
yall = NaN(T+100,n_c);
yall(1:T,:) = yfit_mat(:,1:113);
fall = NaN(T+100,1);
fall(1:T) = ffcst_pct(1:T,3);
fcst = NaN(T+100,5);
fcst(T+1:end,:) = ffcst_pct(T+1:end,:);

% Select Countries for 9 panel plot
%country_list = {'USA';'JPN';'SGP';'ESP';'HUN';'BRA';'CHN';'LBR'};
country_list = {'USA';'JPN';'SGP';'ESP';'HUN';'CHN';'SYR';'LBR'};

figure_4;


% Individual Panel Plots (one for each country)
figure_4_individual_panels;


% Figure 6_b
% Select Countries for 6 panel plot
country_list = {'CAF';'DNK';'GRC';'IND';'KOR';'ECU'};
%country_list = {'CAF';'DNK';'GRC';'IND';'YEM';'ECU'};
figure_6_b;

% Individual Panel Plots (one for each country)
figure_6_b_individual_panels;
